#################################################################################### 
####################################################################################
####################################################################################
# Appendix G: Tables C.18
####################################################################################
####################################################################################
library(here)
library(kableExtra)
library(lmtest)
library(multiwayvcov)
####################################################################################
####################################################################################
## Read in survey data (from: 8_2020_analysis.R)
data = read.csv(here("data/election_panel.csv"))

data$post = 0
data[data$year == 2016 | data$year == 2011,]$post = 1

# Move Mukono into Jinja
data[data$districtname == "Mukono", ]$districtname = "Jinja"

# Create binary indicator for Jinja, because this was the only unit with uneven assignment in the original experiment
data$jinja = 0
data[data$district == 4, "jinja"] = 1

# Create district dummies
for(t in unique(data$districtname)) {
  data[t] = ifelse(data$districtname == t, 1, 0)
}

# Implement SWE of ATE: Demean all covariates, fully interact
cols = c("sex","age","education_level","Arua","Bushenyi","Ibanda","Jinja","Mbale","Mpigi","Pallisa","Mbarara","Shema")
data[, paste0("c_", cols)] = lapply(data[, cols], function(x) (x - mean(x))  )

covariates = c("treatment*c_sex + c_age*treatment + c_education_level*treatment + 
                c_Arua*treatment + c_Bushenyi*treatment + c_Ibanda*treatment + c_Jinja*treatment + c_Mbale*treatment + c_Mpigi*treatment + c_Pallisa*treatment + c_Mbarara*treatment",
               "c_Arua*treatment + c_Bushenyi*treatment + c_Ibanda*treatment + c_Jinja*treatment + c_Mbale*treatment + c_Mpigi*treatment + c_Pallisa*treatment + c_Mbarara*treatment")


# Remove Jinja and Mbarara
datab = data[! data$districtname %in% c("Jinja", "Mbarara", "Mbale", "Bushenyi", "Mpigi"), ]

covariatesb = c("treatment*c_sex + c_age*treatment + c_education_level*treatment + 
                c_Arua*treatment + c_Bushenyi*treatment + c_Ibanda*treatment + c_Pallisa*treatment",
                "c_Arua*treatment + c_Bushenyi*treatment + c_Ibanda*treatment + c_Pallisa*treatment")

# Unrestricted
for (i in names(data[c("election_vote_st","election_attend_st","election_work_st","election_index")])){ # Dependent Variables
  for (j in covariates) {
    model = paste(i,"~","treatment","+","treatment*post","+",j)
    
    # Run each model
    assign(x = paste("m",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = data))
    # Output clustered SEs (county)
    assign(x = paste("c",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = data),
                            cluster.vcov(lm(as.formula(model), data = data), data$village_id)))
  }
}

# Restricted
for (i in names(datab[c("election_vote_st","election_attend_st","election_work_st","election_index")])){ # Dependent Variables
  for (j in covariatesb) {
    model = paste(i,"~","treatment","+","treatment*post","+",j)
    
    # Run each model
    assign(x = paste("bm",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = datab))
    # Output clustered SEs (county)
    assign(x = paste("bc",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = datab),
                            cluster.vcov(lm(as.formula(model), data = datab), datab$village_id)))
  }
}

####################################################################################
####################################################################################
## Table C.18

stargazer(m.election_index.t, bm.election_index.t, m.election_vote_st.t, bm.election_vote_st.t, 
          m.election_attend_st.t, bm.election_attend_st.t, m.election_work_st.t, bm.election_work_st.t,  
          
          se = list(c.election_index.t[,2], bc.election_index.t[,2], c.election_vote_st.t[,2], bc.election_vote_st.t[,2], 
                    c.election_attend_st.t[,2], bc.election_attend_st.t[,2], c.election_work_st.t[,2], bc.election_work_st.t[,2] ),
          
          p = list(c.election_index.t[,4], bc.election_index.t[,4], c.election_vote_st.t[,4], bc.election_vote_st.t[,4], 
                   c.election_attend_st.t[,4], bc.election_attend_st.t[,4], c.election_work_st.t[,4], bc.election_work_st.t[,4] ),
          
          keep = c("treatment$"), 
          
          order = c("$treatment$"), covariate.labels=c("Treatment"),
          
          type = "latex", out = "tables/election_1_c.tex",
          label = "tab:election_1_c",column.sep.width = "1pt", table.placement = "!ht",
          keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
          title = "Effect of LG CHP on Electoral Participation",
          notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", dep.var.caption = "", notes.align = "l", notes.append = F, notes.label = "",
          column.labels = c("Index","Vote","Attend Rally","Work for Party"), column.separate = c(2,2,2,2),
          add.lines = list(c("Restricted", "No", "Yes", "No", "Yes", "No", "Yes", "No", "Yes")))

####################################################################################
####################################################################################
## Table M.49

stargazer(m.election_index.c, bm.election_index.c, m.election_vote_st.c, bm.election_vote_st.c, 
          m.election_attend_st.c, bm.election_attend_st.c, m.election_work_st.c, bm.election_work_st.c,  
          
          se = list(c.election_index.c[,2], bc.election_index.c[,2], c.election_vote_st.c[,2], bc.election_vote_st.c[,2], 
                    c.election_attend_st.c[,2], bc.election_attend_st.c[,2], c.election_work_st.c[,2], bc.election_work_st.c[,2] ),
          
          p = list(c.election_index.c[,4], c.election_index.c[,4], c.election_vote_st.c[,4], c.election_vote_st.c[,4], 
                   c.election_attend_st.c[,4], c.election_attend_st.c[,4], c.election_work_st.c[,4], c.election_work_st.c[,4] ),
          
          keep = c("treatment$"), 
          
          order = c("$treatment$"), covariate.labels=c("Treatment"),
          
          type = "latex", out = "tables/election_1_nc.tex",
          label = "tab:election_1_nc",column.sep.width = "1pt", table.placement = "!ht",
          keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
          title = "Effect of LG CHP on Electoral Participation (No covariates)",
          notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", dep.var.caption = "", notes.align = "l", notes.append = F, notes.label = "",
          column.labels = c("Index","Vote","Attend Rally","Work for Party"), column.separate = c(2,2,2,2),
          add.lines = list(c("Restricted", "No", "Yes", "No", "Yes", "No", "Yes", "No", "Yes")))

rm(list = grep("^(c|m)\\.", ls(), value = T))
rm(list = grep("^(bc|bm)\\.", ls(), value = T))
